N95- 14544 



THERMOCAPILLARY MOTION OF DEFORMABLE DROPS 

Hossein Haj-Hariri and Qingping Shi 
Department of Mechanical, Aerospace, and Nuclear Engineering 
University of Virginia 
Charlottesville, VA 22903, USA 

Ali Borhan 

Department of Chemical Engineering 
The Pennsylvania State University 
University Park, PA 16802, USA 


ABSTRACT 


The thermocapillary motion of initially spherical drops/bubbles driven by a constant temperature gradient 
in an unbounded liquid medium is simulated numerically. Effects of convection of momentum and energy, 
as well as shape deformations, are addressed. The method used is based on interface tracking on a base 
cartesian grid, and uses a smeared color or indicator function for the determination of the surface topology. 
Quad-tree adaptive refinement of the cartesian grid is implemented to enhance the fidelity of the surface 
tracking. It is shown that convection of energy results in a slowing of the drop, as the isotherms get wrapped 
around the front of the drop. Convection of momentum results in a slowdown if the drop has a lower density 
than the exterior phase. Shape deformations resulting from inertial effects affect the migration velocity. The 
physical results obtained are in agreement with the existing literature. Furthermore, remarks are made on 
the sensitivity of the calculated solutions to the smearing of the fluid properties. Analysis and simulations 
show that the migration velocity depends very strongly on the smearing of the interfacial force whereas it is 
rather insensitive to the smearing of other properties, hence the adaptive grid. 


INTRODUCTION 


Recently there has been considerable interest in the development of methods for the numerical simulation of 
flows involving sharp fluid-fluid interfaces. The common feature of the more successful of these methods is 
the replacement of the concentrated interfacial forces by a smeared body force. The discontinuous variations 
of other fluid properties across such interfaces are treated similarly. Therefore the multi-fluid problems 
with concentrated forces reduce to single-fluid ones with smoothly- varying properties, subject to appropriate 
distributions of body forces. The obvious advantage is the simplification of the topology of the problem 
and the requisite numerical grid for the discretization of the equations governing the flow. This nontrivial 
simplification is the direct result of smearing of the properties and the interfacial forces. Some notable 
early contributions are those by Hirt and Nichols [1], Peskin [2], Brackbill et al. [3], followed recently by 
Tryggvason et al. [4], Haj-Hariri et al. [5], and Sheth and Pozrikidis [6]. 

The motion of drops and bubbles in a fluid medium is of fundamental importance in many natural physical 
processes as well as a host of industrial activities. Drops or bubbles floating in a fluid can be moved by a 
large number of ’weak’ forces including gravity, which is normally the dominant mechanism of drop motion 
on earth, and interfacial tension gradients which can play a major role under microgravity conditions where 
sedimentation and gravity-driven convection are largely eliminated. Interfacial tension gradients induced 
at the surface of a drop or bubble through variations of temperature, surfactant concentration, or surface 
charge distribution give rise to unbalanced tangential stresses which result in fluid motion (cf. [7]). 

Most of the work on the motion of drops due to interfacial tension gradients is relatively recent, as sum- 
marized in the review article by Subramanian [8], The majority of the analytical and numerical investigations 
assume a fixed spherical shape for the migrating drop. This is an appropriate assumption only in the limit 
of large interfacial tension with small variations on the surface of the drop. The validity of this assumption 
is highly questionable in the presence of the convective transport of momentum and energy, as well as for 
droplet motion in the vicinity of solid boundaries (e.g. container walls). It is the purpose of this work -and 
others like it- to develop a sufficiently general numerical technique for simulating deformable drops, so that 
effects such as inertia, heat convection, and container walls can be modeled. 
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The thermocapillary migration of droplets was first demonstrated by Young et al. [9], who derived the 
following expression for the migration velocity of a spherical drop of radius a in a constant temperature 
gradient, G, within an unbounded fluid medium: 

Uygb = bG , (1) 

where 6 is a function of the dimensional parameters, whose form is given in the above reference, and is not 
an issue here. Equation (1) is valid for negligible values of the Reynolds and Marangoni numbers, i.e. in the 
absence of convection of momentum and energy. From a practical viewpoint, however, a wide range of values 
of these parameters can be encountered in space applications, depending on the physical properties of the 
fluids involved. Moreover, in practical applications, bubbles and drops interact with each other as well as 
with neighboring surfaces comprising the container walls. This has motivated a series of theoretical analyses 
which have led to considerable progress in predicting the thermocapillary migration velocity of drops at 
nonzero values of Re and Ma, as well as in the vicinity of other surfaces. Szymczyk and Siekmann [10] 
obtained numerical solutions of the equations governing thermocapillary motion of a spherical gas bubble 
for Re < 100 and Ma < 1000, while Balasubramaniam and Lavery [11] performed similar calculations for up 
to Re = 2000 and Ma = 1000. The results of these studies show the scaled thermocapillary velocity to be a 
decreasing function of Marangoni number at fixed Reynolds number. 

Whereas a large range of values of Re and Ma have been explored numerically for the thermocapillary 
motion of a spherical gas bubble, similar studies for the case of liquid drops have been mostly analytical in 
nature, using perturbation techniques which limit the utility of the results to a much smaller range of Re and 
Ma. Moreover, most of the numerical calculations to date have assumed a fixed spherical shape for the drop. 
This limits the usefulness of the large Reynolds number results from these studies since, typically, drops and 
bubbles in common fluids deform substantially when moving at large Reynolds numbers. Some analyses have 
accounted for small deformations of the drop from a spherical shape (e.g. Haj-Hariri et al. [12]). However, 
the results of these asymptotic analyses are limited to a small range of the governing parameters. Since 
the terminal velocity of a drop can be strongly affected by its shape, drop deformations must be taken into 
account. 

In this work the large-Marangoni-number thermocapillary migration of deformable drops and bubbles is 
studied using the method described below. The results are compared with analyses and other comptuations 
where possible. The advantage of the present method is that a deforming drop can be studied just as easily 
as a non-deforming one. Also the case of high-Reynolds-number thermocapillary migration of a deformable 
drop is studied. There is complete agreement with the trends predicted by Haj-Hariri et al. [12]. 


MODELING TECHNIQUE 

Consider the incompressible, axisymmetric thermocapillary migration of a deformable drop in an infinite 
domain in the absence of container walls. The great difficulty in treating this problem lies in the a’priori 
unknown location of the drop surface. The variations in temperature complicate matters further because 
of the intimate coupling between the evolving temperature distribution and the unknown drop shape. We 
account for this explicit coupling by introducing the idea of continuum modeling of the interface between the 
drop and the surrounding phase. The so-called Continuum Surface Force (CSF) model -originally developed 
by Brackbill et al. [3]~ allows the shape of the drop or bubble to be determined as part of the numerical 
solution without any additional effort. 

The basic idea of the CSF technique consists of replacing all surface properties by corresponding volume 
properties defined over a volume obtained by smearing the interface S. This is done using a passive-scalar 
color function C having distinctly different constant values in the two phases. We incorporate temperature- 
induced variation of interfacial tension into the CSF model and modify it to allow for the simultaneous 
solution of the momentum and energy equations. An important observation is that the CSF method fails to 
work for problems involving migration because the tangential component of surface velocities wash the color 
function away from the interface and result in the loss of integrity of the surface after some time. As a remedy, 
the interface C = 0.5, instead of the entire color function field, is advected by the flow. The color function 
is then redistributed based on the new interface location. This is similar to the surface-tracking procedure 
of Unverdi and Trygvasson [4], with the added advantage of grid refinement. Moreover, we suspect that 
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the calcualtion of surface curvatures in three dimensions will be more economical using the color function 
instead of a direct calculation. 

As such, the stress jump, f , across 5, is given by 

where n = — VC/|VC|. The surface force per unit interfacial area is then replaced by its volume-distributed 
’ counterpart, F, where 

F = f jVC| , (3) 

and |VC| plays the role of a smeared delta function. Thus, a single grid can be generated without regard for 
the actual shape of the drop. Following each solution iteration, the color function is propagated passively 
in the computed velocity field, and the deformed shape of the drop is determined. Clearly, the degree of 
smearing is directly related to the resolution of the grid being used. 


GOVERNING EQUATIONS AND SOLUTION STRATEGY 


Including the smeared body force F, expressed by Eq. (3), the dimensionless momentum equation in the 
presence of interfacial tension becomes, 


— 
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( 4 ) 


The solution of the flow problem also requires a knowledge of mass and temperature distributions. The 
mathematical description of the system is then completed by the addition of the dimensionless continuity 
and energy equations : 

V * u — 0 , (5) 

f + " vr =s; v -< 5VT >' < 6 > 

where an overbar denotes the single-phase extension of the physical properties. The variables in the dimen- 
sionless equations above are based on the following characteristic scales, 
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The dimensionless parameters are the Reynolds, capillary, and Marangoni numbers : 
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where cr Q is the interfacial tension at some reference temperature for the drop, and p, p, and a are the density, 
viscosity, and thermal diffusivities of the appropriate phases, respectively. 


The basic numerical algorithm for the solution of Eqs. (4)-(6) is the splitting method consisting of two 
steps. The first step is the calculation of an intermediate velocity field resulting from the convection and 
diffusion of momentum as well as the tangential component of the body force. The second step requires 
the solution of a poisson equation for the pressure to enforce incompressibility. The only explicit step is the 
updating of the interface and the redistribution of the color. 


SENSITIVITY TO SMEARING: ADAPTIVE GRID REFINEMENT 


The method described above is not at its best if it is implemented with inadequate resolution. This was 
shown by Haj-Hariri et al. [13] wherein an analytically tractable model problem was employed to illustrate 
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the effect of the smearing of the fluid properties and the interfacial force on some global variables of the 
problem, in particular the interfacial velocity (analogous to the migration velocity here). This study indicated 
a much stronger dependence of the solution on the smearing of the interfacial forces as opposed to that of 
the other fluid properties. Therefore there is a vital need for enhanced resolution. 

One way to improve the resolution is to use a refined mesh over the interface region. However the use 
of the refined grid in the region away from the interface is unwise because the physical variables have no 
appreciable variations in that region. The adaptive local grid-refinement (quad-tree method) proves ideal 
for this purpose. With the use of local grid-refinement, the computer storage requirements are reduced and 
a desired accuracy is achieved on a near-optimal mesh. An important consideration in adaptive schemes is 
the data structure. The data structure used is a quad-tree structure -“father” cells are refined by division 
into four “son” cells. In this study, whether a cartesian grid cell needs to be refined depends on the gradient 
of the color function C, hence requiring grid refinement only in the vicinity of the interface. Furthermore, 
a smooth transition from large cells to small cells is ensured by prohibiting more than a one-level difference 
between adjacent cells. 


NUMERICAL RESULTS 


Marangoni-number effects 

As a test of the numerical method detailed above, the thermocapillary migration of a spherical drop (bubble) 
with negligible Reynolds and capillary numbers was simulated. The vanishing of these two parameters 
results in a nondeforming drop (bubble) shape. The effect of energy convection (Marangoni number) on the 
migration velocity is presented in Figs. 1 and 2. Figure 1 corresponds to the migration of a gas bubble and 
replicates the work of Balasubramaniam and Lavery [11]. The dimensionless migration velocity for a gas 
bubble in the limit of vanishing Marangoni number is found from Eq. (1) to be 0.5. The increase in the 
convection of energy results in the wrapping of the isotherms about the drop and reduction of the driving 
force and migration velocity (cf. Shi et al. [5]). The migration velocity for this simulation is normalized by its 
value at Ma ~ 0 and presented in Fig. 1. The results are virtually identical to those of Balasubramaniam and 
Lavery [11] for Marangoni numbers less than 100. For higher values of Ma the migration velocity asymptotes 
to 0.255, very close to 0.235, the value predicted analytically by Balasubramaniam [14]. Results from other 
analyses of the same problem, in the absence of any deformations, are presented in table ??. 

If the gas bubble is replaced by a liquid drop (with p = 0.6/w), then for moderate values of Marangoni 
number the effect on the migration velocity is minimal. The migration velocity for the drop is lower than 
that for the gas bubble at all values of Ma. In particular, at Ma — 0, Eq. (1) predicts a value of 0.23 as 
opposed to 0.5 for a gas bubble. 


Reynolds-numb er effects 

In order to allow for deformations, the capillary number was increased to unity. Increasing the Marangoni 
number alone does not seem to lead to any substantial deformations at Re — 0. Hence, we consider instead 
the effects of the convection of momentum by letting Re — 1. The drop phase has density, viscosity, and 
thermal conductivity ratios of 0.6, 0.7, and 0.1, respectively. The migration velocity for such a drop, in the 
limit of vanishing Re and Ca, is given by Eq. (1) as 0.232. The perturbation analysis of Haj-Hariri et al. 
[12] predicts a reduction of the migration velocity, as well as an oblate-spheroidal deformed shape. Both 
of these trends are confirmed by the simulations. The deformed shape is presented in terms of contours of 
the color function in Fig. 3. Figure 4 contains the final adapted grid for this computation. The calculated 
migration velocity of the deformed drop is 0.180. To assess the contribution of the shape deformation to the 
change in the migration velocity, the same simulation was performed while constraining the drop shape to 
remain spherical. The reduction in the migration velocity was now 40% less, at 0.203. Therefore, the role of 
deformations is very important and there is a need for powerful simulation methods that allow the shape of 
the drop to change. 
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CONCLUSIONS 


The model presented above extends previous investigations of the thermocapillary motion of drops by remov- 
ing the assumption of a fixed spherical drop shape and accounting for the convective transport of momentum 
and energy, thereby generating more realistic solutions to this problem. The continuum modeling of free 
interface combined with local grid-refinement provide a powerful tool for studying thermocapillary motions 
of two- and three-dimensional deformable drops. The method eliminates the need for interface tracking and 
reconstruction, and does not impose any restrictions on the number, complexity, or dynamic evolution of 
free interfaces. The deformation of a drop in thermocapillary migration at high values of Re and Ma can 
affect its migration velocity in a manner not predicted by perturbation analysis or fixed-shape simulations. 
In fact, the above results indicate that an increase in the convective transport of momentum and energy can 
lead to drop deformations which, in turn, dramatically affect the thermocapillary migration velocity. 
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Figure 3: Deformed-drop color contours 
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Figure 2: Migration velocity of liquid drop 
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Figure 4: Deformed-drop adapted grid 




